! Bessel J_0(x) function in double precision
!
      function dbesj0(x)
      implicit real*8 (a - h, o - z)
      dimension a(0 : 7), b(0 : 64), c(0 : 69), d(0 : 51)
      parameter (pi4 = 0.78539816339744830962d0)
      data (a(i), i = 0, 7) / 
     &    -0.0000000000023655394d0, 0.0000000004708898680d0, 
     &    -0.0000000678167892231d0, 0.0000067816840038636d0, 
     &    -0.0004340277777716935d0, 0.0156249999999992397d0, 
     &    -0.2499999999999999638d0, 0.9999999999999999997d0 / 
      data (b(i), i = 0, 12) / 
     &    0.0000000000626681117d0, -0.0000000022270614428d0, 
     &    0.0000000662981656302d0, -0.0000016268486502196d0, 
     &    0.0000321978384111685d0, -0.0005005237733315830d0, 
     &    0.0059060313537449816d0, -0.0505265323740109701d0, 
     &    0.2936432097610503985d0, -1.0482565081091638637d0, 
     &    1.9181123286040428113d0, -1.1319199475221700100d0, 
     &    -0.1965480952704682000d0 / 
      data (b(i), i = 13, 25) / 
     &    0.0000000000457457332d0, -0.0000000015814772025d0, 
     &    0.0000000455487446311d0, -0.0000010735201286233d0, 
     &    0.0000202015179970014d0, -0.0002942392368203808d0, 
     &    0.0031801987726150648d0, -0.0239875209742846362d0, 
     &    0.1141447698973777641d0, -0.2766726722823530233d0, 
     &    0.1088620480970941648d0, 0.5136514645381999197d0, 
     &    -0.2100594022073706033d0 / 
      data (b(i), i = 26, 38) / 
     &    0.0000000000331366618d0, -0.0000000011119090229d0, 
     &    0.0000000308823040363d0, -0.0000006956602653104d0, 
     &    0.0000123499947481762d0, -0.0001662951945396180d0, 
     &    0.0016048663165678412d0, -0.0100785479932760966d0, 
     &    0.0328996815223415274d0, -0.0056168761733860688d0, 
     &    -0.2341096400274429386d0, 0.2551729256776404262d0, 
     &    0.2288438186148935667d0 / 
      data (b(i), i = 39, 51) / 
     &    0.0000000000238007203d0, -0.0000000007731046439d0, 
     &    0.0000000206237001152d0, -0.0000004412291442285d0, 
     &    0.0000073107766249655d0, -0.0000891749801028666d0, 
     &    0.0007341654513841350d0, -0.0033303085445352071d0, 
     &    0.0015425853045205717d0, 0.0521100583113136379d0, 
     &    -0.1334447768979217815d0, -0.1401330292364750968d0, 
     &    0.2685616168804818919d0 / 
      data (b(i), i = 52, 64) / 
     &    0.0000000000169355950d0, -0.0000000005308092192d0, 
     &    0.0000000135323005576d0, -0.0000002726650587978d0, 
     &    0.0000041513240141760d0, -0.0000443353052220157d0, 
     &    0.0002815740758993879d0, -0.0004393235121629007d0, 
     &    -0.0067573531105799347d0, 0.0369141914660130814d0, 
     &    0.0081673361942996237d0, -0.2573381285898881860d0, 
     &    0.0459580257102978932d0 / 
      data (c(i), i = 0, 13) / 
     &    -0.00000000003009451757d0, -0.00000000014958003844d0, 
     &    0.00000000506854544776d0, 0.00000001863564222012d0, 
     &    -0.00000060304249068078d0, -0.00000147686259937403d0, 
     &    0.00004714331342682714d0, 0.00006286305481740818d0, 
     &    -0.00214137170594124344d0, -0.00089157336676889788d0, 
     &    0.04508258728666024989d0, -0.00490362805828762224d0, 
     &    -0.27312196367405374426d0, 0.04193925184293450356d0 / 
      data (c(i), i = 14, 27) / 
     &    -0.00000000000712453560d0, -0.00000000041170814825d0, 
     &    0.00000000138012624364d0, 0.00000005704447670683d0, 
     &    -0.00000019026363528842d0, -0.00000533925032409729d0, 
     &    0.00001736064885538091d0, 0.00030692619152608375d0, 
     &    -0.00092598938200644367d0, -0.00917934265960017663d0, 
     &    0.02287952522866389076d0, 0.10545197546252853195d0, 
     &    -0.16126443075752985095d0, -0.19392874768742235538d0 / 
      data (c(i), i = 28, 41) / 
     &    0.00000000002128344556d0, -0.00000000031053910272d0, 
     &    -0.00000000334979293158d0, 0.00000004507232895050d0, 
     &    0.00000036437959146427d0, -0.00000446421436266678d0, 
     &    -0.00002523429344576552d0, 0.00027519882931758163d0, 
     &    0.00097185076358599358d0, -0.00898326746345390692d0, 
     &    -0.01665959196063987584d0, 0.11456933464891967814d0, 
     &    0.07885001422733148815d0, -0.23664819446234712621d0 / 
      data (c(i), i = 42, 55) / 
     &    0.00000000003035295055d0, 0.00000000005486066835d0, 
     &    -0.00000000501026824811d0, -0.00000000501246847860d0, 
     &    0.00000058012340163034d0, 0.00000016788922416169d0, 
     &    -0.00004373270270147275d0, 0.00001183898532719802d0, 
     &    0.00189863342862291449d0, -0.00113759249561636130d0, 
     &    -0.03846797195329871681d0, 0.02389746880951420335d0, 
     &    0.22837862066532347461d0, -0.06765394811166522844d0 / 
      data (c(i), i = 56, 69) / 
     &    0.00000000001279875977d0, 0.00000000035925958103d0, 
     &    -0.00000000228037105967d0, -0.00000004852770517176d0, 
     &    0.00000028696428000189d0, 0.00000440131125178642d0, 
     &    -0.00002366617753349105d0, -0.00024412456252884129d0, 
     &    0.00113028178539430542d0, 0.00708470513919789080d0, 
     &    -0.02526914792327618386d0, -0.08006137953480093426d0, 
     &    0.16548380461475971846d0, 0.14688405470042110229d0 / 
      data (d(i), i = 0, 12) / 
     &    1.059601355592185731d-14, -2.71150591218550377d-13, 
     &    8.6514809056201638d-12, -4.6264028554286627d-10, 
     &    5.0815403835647104d-8, -1.76722552048141208d-5, 
     &    0.16286750396763997378d0, 2.949651820598278873d-13, 
     &    -8.818215611676125741d-12, 3.571119876162253451d-10, 
     &    -2.631924120993717060d-8, 4.709502795656698909d-6, 
     &    -5.208333333333283282d-3 / 
      data (d(i), i = 13, 25) / 
     &    7.18344107717531977d-15, -2.51623725588410308d-13, 
     &    8.6017784918920604d-12, -4.6256876614290359d-10, 
     &    5.0815343220437937d-8, -1.76722551764941970d-5, 
     &    0.16286750396763433767d0, 2.2327570859680094777d-13, 
     &    -8.464594853517051292d-12, 3.563766464349055183d-10, 
     &    -2.631843986737892965d-8, 4.709502342288659410d-6, 
     &    -5.2083333332278466225d-3 / 
      data (d(i), i = 26, 38) / 
     &    5.15413392842889366d-15, -2.27740238380640162d-13, 
     &    8.4827767197609014d-12, -4.6224753682737618d-10, 
     &    5.0814848128929134d-8, -1.76722547638767480d-5, 
     &    0.16286750396748926663d0, 1.7316195320192170887d-13, 
     &    -7.971122772293919646d-12, 3.544039469911895749d-10, 
     &    -2.631443902081701081d-8, 4.709498228695400603d-6, 
     &    -5.2083333315143653610d-3 / 
      data (d(i), i = 39, 51) / 
     &    3.84653681453798517d-15, -2.04464520778789011d-13, 
     &    8.3089298605177838d-12, -4.6155016158412096d-10, 
     &    5.0813263696466650d-8, -1.76722528311426167d-5, 
     &    0.16286750396650065930d0, 1.3797879972460878797d-13, 
     &    -7.448089381011684812d-12, 3.512733797106959780d-10, 
     &    -2.630500895563592722d-8, 4.709483934775839193d-6, 
     &    -5.2083333227940760113d-3 / 
      w = abs(x)
      if (w .lt. 1) then
          t = w * w
          y = ((((((a(0) * t + a(1)) * t + 
     &        a(2)) * t + a(3)) * t + a(4)) * t + 
     &        a(5)) * t + a(6)) * t + a(7)
      else if (w .lt. 8.5d0) then
          t = w * w * 0.0625d0
          k = int(t)
          t = t - (k + 0.5d0)
          k = k * 13
          y = (((((((((((b(k) * t + b(k + 1)) * t + 
     &        b(k + 2)) * t + b(k + 3)) * t + b(k + 4)) * t + 
     &        b(k + 5)) * t + b(k + 6)) * t + b(k + 7)) * t + 
     &        b(k + 8)) * t + b(k + 9)) * t + b(k + 10)) * t + 
     &        b(k + 11)) * t + b(k + 12)
      else if (w .lt. 12.5d0) then
          k = int(w)
          t = w - (k + 0.5d0)
          k = 14 * (k - 8)
          y = ((((((((((((c(k) * t + c(k + 1)) * t + 
     &        c(k + 2)) * t + c(k + 3)) * t + c(k + 4)) * t + 
     &        c(k + 5)) * t + c(k + 6)) * t + c(k + 7)) * t + 
     &        c(k + 8)) * t + c(k + 9)) * t + c(k + 10)) * t + 
     &        c(k + 11)) * t + c(k + 12)) * t + c(k + 13)
      else
          v = 24 / w
          t = v * v
          k = 13 * (int(t))
          y = ((((((d(k) * t + d(k + 1)) * t + 
     &        d(k + 2)) * t + d(k + 3)) * t + d(k + 4)) * t + 
     &        d(k + 5)) * t + d(k + 6)) * sqrt(v)
          theta = (((((d(k + 7) * t + d(k + 8)) * t + 
     &        d(k + 9)) * t + d(k + 10)) * t + d(k + 11)) * t + 
     &        d(k + 12)) * v - pi4
          y = y * cos(w + theta)
      end if
      dbesj0 = y
      end
!
